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ABSTRACT 

Early in the reionization process, the intergalactic medium (IGM) would have been quite inhomogeneous on 
small scales, due to the low Jeans mass in the neutral IGM and the hierarchical growth of structure in a cold 
dark matter Universe. This small-scale structure acted as an important sink during the epoch of reionization, 
impeding the progress of the ionization fronts that swept out from the first sources of ionizing radiation. Here 
we present results of high-resolution cosmological hydrodynamics simulations that resolve the cosmological 
Jeans mass of the neutral IGM in representative volumes several Mpc across. The adiabatic hydrodynamics we 
follow are appropriate in an unheated IGM, before the gas has had a chance to respond to the photoionization 
heating. Our focus is determination of the resolution required in cosmological simulations in order to suffi- 
ciently sample and resolve small-scale structure regulating the opacity of an unheated IGM. We find that a dark 
matter particle mass of m^ m < 50 M and box size of L > 1 Mpc are required. With our converged results we 
show how the mean free path of ionizing radiation and clumping factor of ionized hydrogen depends upon the 
ultraviolet background (UVB) flux and redshift. We find, for example at z = 10, clumping factors typically of 
10 to 20 for an ionization rate of T ~ 0.3-3 x 10~ 12 s _1 , with corresponding mean free paths of ~ 3- 15 Mpc, 
extending previous work on the evolving mean free path to considerably smaller scales and earlier times. 
Subject headings: cosmology: theory — dark ages, reionization, first stars — intergalactic medium 



1. INTRODUCTION 

The fact that the most abundant sources of radiation dur- 
ing reionization were likely to be much less luminous than 
present-day galaxies, combined with their large luminosity 
distances, means that the details of the reionization process 
are beyond most current observational probes. The notable 
exceptions are observations of the polarization of the cosmic 
microwave background (CMB), which i mply an optical depth 
to Thomson scattering of r ~ 0.09 ( jKomatsu et al. 201 l| l, 
and the appearanc e of a Gunn-Pete rson trough in the spectra 
of distant quasars ( Fan et al.][T006). indicating that reioniza- 
tion was largely complete by z ~ 6. Reionization is therefore 
thought to have mainly taken place over the redshift range 
z ~ 6 — 15. Due to the lack of more specific constraints, much 
of our current understanding about the epoch of reionization 
comes from theoretical studies in the context of the ACDM 
cosmology. 

The picture which always emerges is of small-scale gaseous 
structures forming at z > 20, due to the collapse o f dark mat- 
ter halos at the Jeans scale, roughly 10 4 M Q (e.g., [Peebles & 
picke|[T9T8l ICouchman & Rees|[l98l^ |Shapiro et al.||1994{ 
Gnedin & Hui|1998[ l. The gas was just cool enough to fall 
into halos at this mass, leading to strong inhomogenities on a 
scale of tens of comoving parsecs. At the same time, slightly 
more massive halos, with masses on the order of ~ 10 6 M , 
formed enough H2 molecules in their cores to cool efficiently, 
leading to the formation of the first stars in the Universe (e.g., 
Tegmark et aL[[T99"7l |Abel et al.||2"0"02"[ |Bromm et aTJ|2Tj02"' 
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Yoshid a~et al.|2003"[ ). The ionizing radiation from these stars 
is thought to have created substantial, yet short-lived H II re- 
gions, which were sha ped by the surrounding inhomogeneity 
of the gas distributi on ( |Alvarez et al.|2006[|A.bel et aL|2007 
|Yoshida et al.|2007) . 

Eventually, sufficiently large halos formed that triggered 
the formation of the first galaxies , containing tens of thou- 
sands of sta rs Pohnson et al.|2007| |Wise & Abel|2008[ [Grelf| 
|et al.|[2008] >. These nascent dwarf galaxi es would have cre- 
ated longer-lived and isolated H II regions (Wise & Cen 2009, 
[Wise et aL]|2012| ). It is unclear how these galaxies evolved 
into the much more luminous o nes that have been observed at 
redshifts as high as z ~ 8 (e.g., |Bouwens et aT7||2010 1. Nev- 
ertheless, it is widely believed that as the first galaxies grew 
and merged, their collective radiative output created a large 
and complex patchwork of ionized bubbles, with characteris- 
tic sizes on the order of tens to hundreds of comoving Mpc 



(e.g.J Barkana & Loeb 2004; Furla netto et al.|2004[|lliev et al. 
|2006f . During this time, dense systems i n the IGM likely 



impeded the progress of ionization fronts (Barkana & Loeb 
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1999, Haiman et al. 2001; Shapi ro et al]|2004| |Iliev et~aL 
2005 ; Ciardi et al. 2006). At the end of reionization the so- 
called "Lyman-limit" systems, dense clouds of gas optically- 
thick to ionizing radiation observed in the spectra of quasars 



at z < 6 (e.g., |Storrie-Lombardi et al.] [1994; Prochask a et al. 
2009 ), dominated the overall opacity of the IGM to ionizing 
radiation. These systems crucially influenced t he percola- 
tion phase of reionization (|Gnedin & Fan 2006, Cho udhury] 
|et al.|2009[|Alvarez~& Abel 2 0l2| l, which in turn determined 
the evolution and structure of the ionizing background (e.g., 
I Haardt & Madau[[T996[ |Bolton & Haehnelt|[2007) [McQuinn 
|etal.|2011) .~ 

Thus, the progress of reionization depended not only on the 
properties of the sources of ionizing radiation, but also on the 
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sinks. Theoretical models of reionization must describe not 
just the spectral energy distribution, abundance, and cluster- 
ing of early sources of ionizing radiation, but also the inho- 
mogeneity of the intergalactic medium (IGM) in the space 
between the sources. It is this latter description that is the 
goal of the present work. 

Early descriptions of reionization took into account inho- 
mogeneities in the IGM through a "clumping factor", c;, by 
which the recombination rate is boosted relative to the homo- 
geneous case. This allows one to write a global ionization 
rate, equal to the ionizing photon emmissivity minus the re- 
combination rate of a clumpy IGM, and thereby determine the 
reionization history for a given ionizing source population. 
Shapiro & Giroux ( 1987) used such a model to show that the 



TABLE 1 
Simulation Parameters 



observed population of QSOs were insufficient to have reion- 
ized the Universe by z ~ 5. Their assumption of c\ ~ 1 would 
have been conservative, in that that additional recombinations 
would have made it even more difficult for quasars alone to 
reionize the Universe. 

In addition to being useful in modeling the reionization his- 
tory, the clumping factor is also important in estimating the 
necessary number of ionizing photons per baryon to maintain 
an ionized Universe. The necessary and sufficient condition 
for maintaining an ionized Universe is that the ionizing pho- 
ton emissivity should be greater than or equa l to the recom- 
bination rate of the IGM. |Madau et"aL] ( |1999| ) used this fact 
to derive a critical star formation rate, above which the rate 
of ionizing photons is enough to maintain the Universe in an 
io nized state. 

[Gnedin & Ostriker (1997]) used hydrodynamic simulations 
with a treatment of photoionization in the "local optical 
depth" approximation to determine the clumping factor of the 
ionized component of the IGM, finding a value of c; ~ 30 at 
2 = 6. They also pointed out that the actual clumping factor of 
the IGM would have been larger due to stru cture on s maller 
scales than they resolved. More recently, Miralda-Escude 
et al. (2000} built a semi-analytical model for the reioniza- 



tion of an inhomogenous IGM, in which the underlying gas 
density distribution was determined by numerical simulations. 
They argued that in addition to specifying the clumping factor 
of the ionized medium, it is also necessary to describe the dis- 
tribution of high-density gas clouds that are able to self-shield 
agains t ionizing r adiation. 

McQuinn et al.|(|2011[) follow ed a similar approach to that 
of Miralda-Escude et al. (2000) to explain the evolution of 
the ionizing background radiation at redshifts less than z ~ 6, 
using more realistic numerical simulations which were post- 
processed with radiative transfer. These works were focused 
on the large scales relevant in the post-reionization IGM, after 
photoionization heating has "ironed out" the dumpiness of 
the IGM on the smallest scales. Our work here is focused on 
the higher redshifts and smaller scales that were most relevant 
early in reionization. 

We sought convergence to the measures of inhomogeneity 
of the unheated IGM during the epoch of reionization, like 
the mean free path, A, clumping factor, q, and density thresh- 
old above which gas is self-shielded, A cl ; t , by spanning the 
parametre space of redshift and ionizing background intensity, 
j u . To do this, we post-process cosmological adiabatic hy- 
drodynamics simulations with radiative transfer calculations 
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along different lines of sight through the simulated volume. 
Radiative feedback raises the Jeans mass of the IGM, thereby 
increasing the scale of inhomogeneities. Therefore, the reso- 
lutions we find in our adiabatic simulations that are necessary 
to resolve structure in the unheated IGM are also sufficient to 
model radiative feedback at all times. 

The outline of the paper is as follows. Details of the sim- 
ulation setup and radiative transfer are described in §2. In §3 
we present our numerical results, followed by §4, where we 
present the results of our convergence tests. §5 concludes with 
a discussion of our main results. 

2. NUMERICAL APPROACH 

Here we describe our numerical approach, in which we 
perform a suite of cosmological adiabatic hydrodynamics 
simulations using the publicly available SPH code Gadget- 
2 (Springel 2005 ). We then postprocess each simulation with 
multifrequency radiative transfer of hydrogen ionizing radia- 
tion, assuming photoionization equilibrium, to determine the 
dependence of basic quantities, like the ionizing photon mean 
free path and clumping factor, on redshift and intensity of the 
background radiation field. 

2.1. Cosmological Hydrodynamic Simulations 

The cosmological simulations are parameterized by box 
size, L, and total number of dark matter and gas particles, N. 
Table [T] summarizes these parameters for our suite of simula- 
tions and lists their corresponding dark matter particle masses, 
radm, along with the comoving gravitational softening length, 
r SO ft. The simulations were evolved from redshift z = 200 to 
z = 6, except for simulations CI through C4 which, due to 
computational limitations, were terminated early at z = 10. 
Initial conditions were generated separately for dark matter 
and baryons using transfer functions computed by CAMB for 
each component, with the same random phases. Through- 
out our work we assume the set of cosmological parameters 
(fi DM , Q b , fl A , h) = (0.228, 0.042, 0.73, 0.72). 

A quantitative test of the simulated structure formation is 
to identify dark matter halos to construct mass functions, 
dn/dM, which can be compared to analytic models. FigurefT] 
shows the mass functions obtained from a friends-of-friends 
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FIG. 1. — Halo mass functions from six simulations are compared to their 
Warren et al, counterparts using the variance £^ e g — cr — ^~(.^bo\) &t redshifts 
10 and 20. The top panel plots simulations Al, B2, and C3 each containing 
dark matter masses of mj m = 31 Mq while the lower panel plots simulations 
A4, B5, and C6 with = 1.6 X 10 4 Mq. In each case, points denote 
halo mass functions obtained from the simulations while the lines trace the 
corresponding Warren et al. curves. 



(FOF) halo identification scheme with linking length of 0.2 
mean interparticle spacings, at redshifts z = 10 and 20 for two 
groups of simulations sharing common mass resolutions of 
'«dm = 31 Mq and 1.6 x 10 4 M . A co mmon fitting function 
to compare to is the | Warren et a l. (2006) mass function. When 
doing so, however, it is important to note that this model as- 
sumes a Universe with infinite spatial extent; something that 
cannot be achieved using numerical simulations. It is there- 
fore useful to compute Warren et al. mass functions using 
a modified variance of the form o% s = a 2 - a 2 (M\, ox ), where 
Mbox is the total mass contained within the simulated volume. 
This has the effect of removing contributions from mass fluc- 
tuations outside on scales larger than that of the simulated vol- 
ume. With this correction we find that the Warren et al. mass 
function is generally well matched by the numerical simula- 
tions. 

There is one important feature worth noting in Figure [T] 
For fixed mass resolution, simulations with a larger number 
of particles (i.e., larger volumes) tend to trace the analytic 
curves more closely. This is most noticeable in the top panel 
for z = 20. We can attribute this to the fact that at fixed res- 
olution, simulations with larger volumes will contain a more 
statistically representative collection of halos. In §4 we will 
show how sample variance in small boxes has important con- 
sequences for numerical convergence. Even though a simula- 
tion may have a sufficient mass resolution to resolve low-mass 
halos within the IGM, its volume may be so small that sample 
variance becomes important when interpreting its results. 

It is well known that the inhomogeneous nature of the IGM 
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FIG. 2.— Gas PDFs for our fiducial simulation B2 (512 3 , 0.5 Mpc) at dif- 
ferent redshifts, as labelled beside the curves. The vertical dot-dashed line 
delineates the overdensity A v i r = 18-7T 2 of a virialized isothermal sphere. If 
collapsed structures have a density profile of the form p oc r , for A v jj < 
p/p < A max , then P(A) oc A -2 5 for A vir < A < A max . We plot the PDF 
multiplied by A 2 s so that the curve should approach a constant for A > A v j, 
if gas is collapsed within these structures. 

plays an important role in the progression of the reionization 



epoch. This was emphasized by Miralda-Escude et al. ( 2000 ) 
who presented an evolutionary model of reionization based 
on the gas density distribution observed in numerical simula- 
tions. It is therefore useful to examine the density distribution 
of baryons within cosmological simulations, through the use 
of the probability density function (PDF), P(A), defined to 
be the normalized distribution of gas in terms of overdensity 
A = P /p. 

In Figure [2] we plot volume-weighted gas PDFs from our 
fiducial simulation B2 which contains 2 x 512 3 dark matter 
plus gas particles in a box of comoving length 0.5 Mpc. More 
precisely, we plot A 2 5 P(A) which is expected to approach a 
constant at A > A v ; r = 187r 2 if gas at those densities is col- 
lapsed within halos described by a density profile p oc r~ 2 . As 
time evolves, the fraction of gas collapsed within halos in- 
creases, though only for z < 10 does the high-density tail of 
the PDF appear to approach a constant value. Note that the 
PDFs shown here are most appropriate to the unheated IGM, 
and after substantial reheating, we expect the distribution to 
evolve in such a way as to decrease the amplitude of P(A) at 
large values of A. 

2.2. Post-processed Ionization Calculation 

In order to simulate the effects of self-shielding by ab- 
sorption systems, we postprocess the SPH density field with 
a multifrequency radiative transfer algorithm. This involves 
tracing the attenuation of the ionizing radiation along differ- 
ent lines of sight throughout the volume while assuming pho- 
toionization equilibrium. 
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2.2. 1 . Ultraviolet Background Spectrum 

We consider a background ionizing intensity I v , so that the 
flux of photons capable of ionizing hydrogen is 



dfl 



hv 



dv, 



(1) 



where hv m = 13.6 eV is the photon energy at the Lyman edge. 
The upper limit in the integral corresponds to the ionizing 
threshold for fully ionizing helium - we assume helium is 
singly ionized along with hydrogen, and therefore only con- 
sider photons below the He II Lyman edge. 
We adopt a power-law UVB spectrum, 



h=Io( — 



(2) 



where vhi < v < Avhi and I a is the intensity at the Lyman 
edge. In our analysis we have sampled a region of parame- 
ter space for which 1 < a < 3. Our fiducial value of a = 2 
is chosen to be consistent with the spectral index we would 
expect for an ionizing backgr ound produced from a m ixture 
of galaxies and quasars (e.g., Bolton & Haehnelt 2007J. Our 
results exhibit only a minor dependence on spectral index in 



this range, as also found by |McQuinn et al.] l |2011) . For this 
reason we henceforth refer only to our fiducial case of a = 2. 

The intensity is often expressed in terms of the quan- 
tity 7_2i, defined to be the isotropic equivalent of I„, 
(J 7 ^^)/(47r), in units of 10" 21 erg cnrV'Hz^ster" 1 . For 
the form expressed in equation |2]), we can integrate equation 
(fill to relate J-21 to the flux of ionizing photons. For a = 2, we 
obtain: 

7-21= 0.09 f F - 2 . l )- (3) 
\ 1(P cm z s 1 J 

Another useful quantity to describe the UVB is T-n, defined 
to be the ionization rate per atom, in units of 10~ 12 s" 1 : 



r_ 12 = 0.3 



10 5 cm 2 s 



(4) 



Note that this refers to the ionization rate corresponding to a 
given background, and not the mean ionization rate per atom 
along our rays, which is lower due to attenuation and includes 
the neutral component of the IGM. 

We calculate the ionization state of the volume for a broad 
range of background flux with the fiducial value of F = 
10 5 cm" 2 s _I taken to be consistent with the value of T_i2 ~ 0.3 
inferred from the optic al depth of the Lyq fo rest seen in 
quasar spectra (e.g., Bolton & Haehnelt 2007). Due to its 
common usage in the literature, we will report our results in 
terms of T_i2, though it should be remembered that its con- 
version to flux simply follows equation 

2.2.2. Ray-tracing 

In our ray-based approach, the UVB has a plane-parallel 
direction dependence, so that I u = F„<5(n), where h is the di- 
rection of propagation of the radiation and F v is the spectral 
flux density. This is appropriate especially in the beginning 
stages of reionization, where a given patch of the IGM is ini- 
tially exposed to a one-sided flux from the downstream di- 



rection of the ionization front. In addition, because recom- 
binations to the ground state are quickly canceled by subse- 
quent photoionizations, we treat the rays independently and 
use the "case-B" recombination coefficient. Finally, because 
the equilibration time is very short compared to the Hubble 
time, we use photoionization equilibrium, which allows us 
to calculate the ionization state and attenuation of the back- 
ground self-consistently by sequentially iterating along the 
ray in the direction ft. 

To obtain an unbiased sample of the gas density field and 
minimize noise, the rays are assigned starting points uni- 
formly distributed in a plane with orientations perpendicu- 
lar to the plane. We use three orthogonal planes in order to 
sample different directions. Each ray segment corresponds 
to a cubic volume element, within which the mean density is 
obtained from the SPH particl e data by the m ass-conserving 
spline interpolation outlined in |Alvarez et al.| ( [2006| ). The ray 
segments have lengths given by L/N my , where L is the box 
size, so that the number of rays is proportional to 7V 2 ay , while 
the number of segments along a given ray is proportional to 
N my . We check for convergence in our radiative transfer cal- 
culations by interpolating to a variety of values for 7V ray . From 
this we find that it is necessary to interpolate to N iay = 1024 
for the 256 3 and 512 3 particle simulations and to iV ra y = 2048 
for the 1024 3 particle simulations. 

2.2.3. Equilibirum Radiative Transfer 
The equation of radiative transfer for l v is 
dl v 



ds 



(5) 



where s is the proper distance and «hi is the proper number 
density of neutral hydrogen. Here we are concerned with the 
transfer of ionizing radiation through a patch of IGM in which 
there are no sources, and therefore set e v = 0. 

The ionization rate within a given ray segment is related to 
the background intensity and neutral hydrogen column den- 
sity, Nhi, integrated along a given ray by using the solution of 
equation (T5]): 

l-Av m 1 



r = 47T 



hv 



-Nhi a 1 



dv, 



(6) 



where u v is the absorption cross section, with mean value 



(■AVHI /„ 
/ n - 

— ■'vhi 



\ Avm hdv 

J Vhi hv 



: 2.84 x 10" 18 cm 2 



(7) 



The flux of ionizing photons is diminished along the ray 
through absorption by intervening neutral hydrogen, and the 
spectrum steepens as softer photons are preferentially ab- 
sorbed. 

To determine the opacity along the ray self-consistently, we 
iterate along the ray, using the total H I column density from 
the previous ray segments to calculate the photoionization rate 
at the current segment using equation (|6}. This is then used to 
determine the neutral hydrogen density in the current segment 
under the assumption of photoionization equilibrium: 



Trim = OL B n HU n e 



(8) 
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where rini, rinu and n e are the number densities of neutral hy- 
drogen, ionized hydrogen, and electrons within that segment. 
The resulting value of H I density is then used to update the 
total column density, and the procedure is repeated until the 
end of the ray is reached. 

Equation dSl assumes a uniform radiation field within each 
ray segment. This assumption breaks down if the segment be- 
comes sufficiently optically thick that Y changes significantly 
across it. To address this issue, individual ray segments are 
split into plane-parallel subsegments in the direction h with 
widths chosen such that the flux passing through each subseg- 
ment is attenuated by no more than 2% of its initial value. 
Photoionization equilibrium is applied in sequence to each 
subsegment and global quantities pertaining to the segment 
as a whole are computed as volume averages over each sub- 
segment. 

We assume that all free electrons within the volume come 
from hydrogen and consider a uniform gas temperature of 
7g as = 10 4 K so that a B = 2.6 x 10" 13 cmV. Including he- 
lium in our calculations would lead to small corrections in the 
hardening of the radiation at high optical depths, due to the 
slightly different frequency dependence of the He I absorption 
cross-section relative to that of H I. Given the insensitivity of 
our results to varying the spectral slope a, inclusion of he- 
lium radiative transfer would not improve the accuracy of our 
results, while needlessly complicating their interpretation, so 
we neglect it. 

2.2.4. Optimal Ray Length and the Mean Free Path 

The opacity of the IGM can be written in terms of the mean 
free path, with the equation of transfer for the flux of ionizing 
photons given by 

as A 

where F is the total flux of ionizing photons in units of 
cm~ 2 s _I , and the factor l+z accounts for the fact that we de- 
fine the mean free path to be in comoving units, and is the def- 
inition we use throughout this paper. We calculate the mean 
free path along a given ray as the solution to equation |9]i: 

A = ~ . * /F . , (10) 

where s is now the comoving length of the ray, F is the inci- 
dent flux at the start of the ray, and F out is the attenuated flux 
leaving the last ray segment. The total mean free path is de- 
termined by first averaging F out over all rays, then applying 
equation ( 10 1. 

Naively, we may choose to set s = L so that each ray sam- 
ples the entire length of the box. However, we must be careful 
since the use of equation ( fT0) i is not physically meaningful in 
the optically thick limit where F out can tend to 0. In other 
words, we want to determine the opacity of the IGM due to 
small-scale structure at a fixed background flux, but includ- 
ing the cumulative effect of this shielding over distances ap- 
proaching the mean free path would correspond to a lower 
flux than what we assume. 

An easy way to avoid this problem is to send photons along 
shorter rays. Of course, this has the disadvantage of sam- 
pling smaller portions of the IGM, possibly missing individ- 
ual self-shielded structures. It is thus optimal to choose a ray 
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FIG. 3. — Here we demonstrate our procedure for choosing the optimal 
ray length that adequately samples the IGM while remaining optically thin. 
The data pertains to simulation B6 (512 3 , 8 Mpc) taken at redshift z = 10. 
The nearly horizontal lines show A as a function of ray length s for different 
photoionization rates where the labels denote T_i2 in units of 10~ 12 s _I . For 
each of these lines, the optimal ray length is chosen as the largest value of s 
for which i < X(s)/5 and is denoted by a black cross. For comparison, the 
diagonal line traces out A = i so that portions rightward of this curve belong 
to the optically thick regime where the mean free path changes significantly. 

length such that the rays on average remain optically thin, 
while still sampling a sufficiently long distance to take into 
account the self-shielding of individual dense structures. We 
achieve this by first calculating A as a function of s, and then 
choose the optimal ray length to be the largest value of s for 
which s < \(s)/5. A lower cutoff of s > L/32 is also applied. 
If this condition cannot be satisfied we flag the given region 
of parameter space and omit its inclusion in our analysis. For 
ray lengths s < L, the usage of the box is maximized by reset- 
ting the intensity along each ray after a distance of s, until the 
ray has traversed a distance L. 

Figure [3] demonstrates our procedure of selecting ray 
lengths atclifferent fluxes for simulation B6 at z = 10. The 
mean free path converges in the optically thin limit where 
A > s but begins to deviate strongly during the transition to 
the optically thick transition when A approaches s. It is clear 
from the plot that an erroneous value for A would be obtained 
for an improper choice of s. The convergence of A in the 
optically thin limit shows that our choice of picking s to be 
bounded by A/5 is robust in changing this fraction by a factor 
of a few. 

3. SIMULATION RESULTS 

We first describe the dependence of clumping factor, q, of 
ionized gas on redshift and photoionization rate. Next, we 
use the gas PDF matched to the clumping factors we have 
obtained, to define a critical overdensity, A cr i t , above which 
gas remains self-shielded and neutral. The opacity of the IGM 
to ionizing radiation, expressed in terms of the mean free path, 
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A, is described next. We first discuss its overall properties 
and then show how it can be used to relate the emissivity of 
ionizing sources to the photoionizing background that they 
produce. Finally, we compare the clumping factors and mean 
free paths obtained here to those which would be expected for 
an optically thin model of the IGM. 

The main goal of this work is to assess the small-scale con- 
vergence of numerical quantities during the initial phase of 
reionization. This is presented in §4 where we show simu- 
lation C3 (1024 3 , 1 Mpc) to be our "converged" simulation. 
However, since C3 was only evolved to z = 10, we show here 
results from our fiducial simulation B2 (512 3 , 0.5 Mpc) in or- 
der to present results down to z = 6. Table [2] summarizes the 
clumping factors, critical overdensities, and mean free paths 
at select redshifts and photoionization rates for simulation B2. 
These values are within 6% of those from the converged sim- 
ulation C3 for all z > 10. 

3.1. Clumping Factor 

Studies of reionization typically make use of the clumping 
factor of ionized gas, defined as 



TABLE 2 

Opacity of the Unheated IGM at Select Values 



Ci 



(nl)/(n e 



(11) 



where n e is the number density of free elections and angled 
brackets denote volume averages over space. The clumping 
factor describes the enhancement of the recombination rate 
relative to a uniform gas distribution, and is therefore crucial 
in understanding the role of inhomogeneities in the ionizing 
photon budget during and after reionization. 

In Figure |4] we plot the clumping factor obtained from the 
radiative transfer calculations performed on our fiducial sim- 
ulation B2. Some general trends of the clumping factor are 
shown by comparing the two panels of this figure. In the first 
place, at fixed redshift, c/ increases with the strength of the 
ionizing background as the flux of ionizing photons are able 
to penetrate further into thick gas clouds, exposing their dense 
interiors where the recombination rate is greatest. Eventually, 
the clumping factor plateaus at large flux where the ionization 
state of the box has saturated and all the gas is ionized. At this 
point ci tends to the total clumping factor of gas in the box as 
n e approaches n. Furthermore, at fixed photoionization rate 
we observe c/ to increase with decreasing redshift which is 
consistent with ongoing structure formation within the IGM. 

Historically, clumping factors of q ~ 30 at z ~ 6 hav e 
been found to be appropriate (e.g., |Gnedin & Ostriker|1997| l, 
though more recently there has been a growing trend towards 
values an order of magnitude smaller. It thus appears contrary 
to historical development that we reproduce c\ ~ 30 at z = 6 
with our fiducial case of T_i2 = 0.3, and find even l arger values 
with increased flux. However, as explained by Pawlik et al.| 



(2009), the passage of an ionization front through the IGM 
will photoevaporate the smallest halos in the box and conse- 
quently suppress the evolution of the clumping factor at small 
scales as the gas is dispersed back into the diffuse IGM. Since 
we do not include such hydrodynamic feedback processes in 
our analysis, the values reported here cannot be used in refer- 
ence to an IGM that has been heated through photoionization. 
Nevertheless, our values are perfectly applicable to the early 
stages of reionization, before the gas has had time to respond 
to the ionizing radiation field. Moreover, as discussed in §5, 
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previous simulations have underestimated the clumping fac- 
tor by a factor of ~ 3 during this period, and may therefore be 
underestimating its subsequent evolution and the impact that 
unresolved small-scale structure had in regulating the early 
stages of reionization. 

3.2. Critical Overdensity for Ionization 

Since the clumping factor describes the distribution of ion- 
ized gas within the volume, it is in principle derivable from 
knowledge of the gas PDF and details of the photoionizing 
radiation field. In a simplified description, we assume that 
all gas within the box with overdensity A < A cl -i t is ionized, 
while the rest is neutral. This is obviously an idealized de- 
scription of reality where a gradual transition between ionized 
and neutral regions will necessarily occur. Any departures 
from the simplified model reflect variations in the local ion- 
izing background and degree of self-shielding and shadowing 
within the inhomogeneous IGM (Miralda-Escude et al. 2000). 

In the simplified model the clumping factor of ionized gas is 
related to the total gas PDF through the following expression: 



ci ■ 



/ () A -A 2 P(A)JA 
(f A ""AP(A) dA 



(12) 



where A cr j t is interpreted as the critical overdensity above 
which self-shielding prevents the gas from becoming ionized. 
It is often useful to assume the form of equation ( p"2} taken 
with some nominal choice for A cr j t in order to compute q 
from a given gas PDF. For example, Chiu et aL| ( |2003| l con- 
sider a model where all gas within collapsed halos is self- 
shielded while all remaining gas is subject to ionization from 
an UVB. In this case, A cr i t = 6tt 2 corresponding to the over- 
density at the virial radius of an isothermal sphere with a mean 
overdensity of A v ; r = 187r 2 . 

Since we compute the clumping factor directly from our ra- 
diative transfer calculations, we take the opposite approach, 
inverting equation ( 12 1 in order to compute A cl i t from knowl- 
edge of P(A) and c/. In doing so, we observe the expected 
trend that A cr i t increases when the photoionization rate is in- 
creased, making the medium more susceptible to ionization. 
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FIG. 4. — (left) Clumping factor versus photoionization rate T_i2 in units of 10 12 s 1 for different redshifts, as labelled beside the curves, (right) Clumping 
factor versus redshift for different photoionization rates with labels denoting the value of r_i2. In both cases the data pertains to our fiducial simulation B2 (512 3 , 
0.5 Mpc) and only regions of parameter space satisfying the ray length criterion described in §2.2.4 are shown. 



In fact, we find the rough proportional ity A c r j t oc r_ 12 which, 
from equation (4) of McQuin n et al.| ( |201 is expected for 
a PDF satisfying P(A) oc A~". We showed in Figure [2lthat 
our PDFs satisfy this power-law at z < 10 for A > A v ; r . This 
is consistent with the model where gas at these densities is 
collapsed in isothermal spheres. Around our fiducial value of 
r_i2 = 0.3, we further find that A cr i t is roughly proportional to 
(1 + z)~ 3 , indicating that the critical proper hydrogen number 
density, n cr j t , is rather insensitive to redshift. A good value to 

take is w C nt ~ 0.1 cm" 3 T^ 3 ,. 

The validity of the idealized model where all gas with over- 
density A < A cr ; t is ionized is tested in Figure [6] Here we 
plot the total gas PDF along with the ionized and neutral 
PDFs obtained from our radiative transfer calculation using 
our fiducial parameters T_i2 = 0.3 and z = 10 for simulation 
B2. The vertical dot-dashed line shows the corresponding 
value of A c rit - its role in delineating the neutral and ionized 
portions of the gas is clearly visible. As anticipated, the tran- 
sition between ionized and neutral regions is not sharp, but 
rather gradual as a consequence of the spatially varying ion- 
izing background and self-shielding due to dense gas pock- 
ets. Nevertheless, our findings indicate that the approximation 
that A cl it represents the critical overdensity above which self- 
shielding maintains the neutral state of the IGM is generally a 
good one. 

3.3. Mean Free Path 

We quantify the opacity of the IGM to the exposed UVB 
through the use of the mean free path of ionizing radiation. 
Conceptually, one can consider the mean free path to be af- 
fected by two components: a diffuse gaseous phase that per- 



vades the IGM and thick gas clouds embedded within col- 
lapsed dark matter halos. The latter make up a significant frac- 
tion of absorption systems that have neutral hydrogen column 
densities Nhi <; l/cr ~ 10 17 cm 2 allowing them to self-shield 
against ionizing radiation. It is within these optically-thick 
structures where the global recombination rate of the IGM is 
dominated and the majo rity of ionizing photons are absorbed 
(Miralda-Escude et al. 2000). As a result, they can signifi- 
cantly impede the progress of reionization. 

The mean free path obtained from our radiative transfer cal- 
culations is computed through the use of equation ( [T0| which 
naturally encompasses both the clumpy IGM and nalo com- 
ponents. In Figure [5] we plot the mean free path as a function 
of photoionization rate for fixed redshift and also as a func- 
tion of redshift for fixed photoionization rate. Some general 
trends are immediately clear in this plot. In the first place, at 
fixed redshift we see that the mean free path increases with the 
strength of the ionizing background. A stronger flux of ioniz- 
ing photons will naturally penetrate further through a diffuse 
IGM and overcome thicker self-shielding structures, consis- 
tent with the previous observation that A cr it oc T^j 3 ,. In addi- 
tion, when the ionizing background is held constant, the mean 
free path is found to increase with decreasing redshift. 

It is important to note that there are two competing fac- 
tors affecting the redshift evolution of A. On the one hand, 
the expansion of the Universe continually dilutes the density 
of hydrogen, hence favouring a strong increase in A with z. 
On the other hand, increased structure formation at low red- 
shift enhances the distribution of Lyman-limit systems that 
strongly inhibit the distance an ionizing photon can propagate 
through the IGM before being absorbed. In the right panel 
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FIG. 5. — (left) The mean fee path versus photoionization rate T_i2 in units of 10 12 s 1 for different redshifts, as labelled, (right) Mean free path versus redshift 
for different photoionization rates with labels indicating the value of T_i2. In both cases the data pertains to our fiducial simulation B2 (512 3 , 0.5 Mpc) and 
only regions of parameter space satisfying the ray length criterion described in §2.2.4 are shown. The dott ed b lack curves in the right panel show the mean free 
path expected for an optically thin, completely ionized, and homogeneous IGM as expressed in equation j!7fr . From bottom to top in the plot, the dotted lines 
take r_!2 = 0.03, 0.3, and 3, and are each calculated using x = Ci = 1. We expect the dotted lines to converge with our results at high redshift when the medium 
approaches homogeneity. At low redshift we observe a large suppression in the calculated mean free path that results from increased structure formation within 
the inhomogeneous IGM. 



of Figure [5] we compare A obtained here to equation ( 17 i for 
the same set of photoionization rates. Taking x = c; = fin this 
equation yields the mean free path we would obtain in an op- 
tically thin, homogeneous, and completely ionized medium. 
In such a model the mean free path evolves rapidly with red- 
shift as A oc (1 +z)~ 5 . Instead, we observe a strong suppres- 
sion in mean free path at low redshift compared to equation 
( fT7| . This highlights the important contribution from inho- 
mogeneities in the IGM. 

3.4. Relationship Between Emissivity and Photoionization 

Rate 



In the context of reionization it is desirable to know the 
ionizing background produced by some population of sources 
with known emissivity. This relationship can be found by 
solving 



r=« ion (i+z)'A(r,z)<7, 



(13) 



where n\ on is the comoving ionizing emissivity and A(F,z) is 
the comoving mean free path that depends on both the ion- 
izing background and redshift. In Figure [71 w e show the de- 
pendence of r on «i on by solving equation ("13 i with the mean 
free paths taken from our radiative transfer calculations. We 
find that T exhibits a rather steep dependence on emissivity 
and appears to diverge at large values of «j on . This behaviour 
is attributed to the fact that not only are there more ionizing 
photons as the emissivity rises, but also their ability to pene- 
trate further through the IGM increases. 
We can relate this behaviour back to the dependency of A on 



r. For instance, suppose we have the simple relation A oc T 13 
at some redshift. Then from equation jl3| we will have that 
F oc ri7 Qn where 7 = (1 — p)" 1 . In Table Bjwe list the values of 
7 obtained by fitting a power-law to our fiducial mean free 
paths within the range 0.1 < r_n < 1 at different redshifts. 
This flux range is considered to emphasize the the relationship 
between T and n; on around our fiducial value of T_i2 = 0.3. We 
find that the relationship between T and «; on strengthens as the 
redshift increases - 7 varies from 2.6 to 14 between redshifts 
6 and 20 respectively. This occurs because the slope f3 rises as 
the IGM becomes more uniform, approaching a limiting value 
of unity for a completely homogeneous Universe with q = 1 
in equation ( flT} - a manifestation of "Olber's Paradox". 

From this trend we can deduce that decreasing the simula- 
tion resolution should steepen the curves in Figure [7J as the 
density distribution becomes more homogenous. Indeed, this 
relation is observed between our suite of simulations where 
we find 7 = 3.2 at z = 6 for our worst-resolved simulatio n A6 



Mc- 



The 



(256 3 , 8 Mpc ), compared to 7 = 2.6 for simulation B2. 
Qui nn et al.| ( |201 1) report the value of 7 ~ 4 at z = 6. 
discrepancy with our result likely arises from a combination 
of our increased resolution and our omission of photoheating 
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FIG. 6. — The bottom panel compares the total gas PDF of our fiducial sim- 
ulation B2 (512 3 , 0.5 Mpc) to the PDFs of neutral and ionized gas within 
the box after radiative transfer is applied. Ray segments are labelled as ion- 
ized if they have an ionized fraction greater than 0.5 and are labelled neutral 
otherwise. The data corresponds to a snapshot at z = 10 with an ionizing 
background of T_i2 = 0.3. The vertical dot-dashed line denotes A cr [ t deter- 
mined by comparing c/ obtained fro m the radiative transfer calculation to the 
total gas PDF through equation \ \2) . The top panel shows the corresponding 
volume fraction of neutral and ionized gas as a function of A. 




10 61 10 s * 
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FIG. 7. — The photoionization rate, expressed in terms of T_i2, as a func- 
tion of the comoving ionizing emi ssiv ity at different redshifts. These curves 
are evaluated by solving equation l |13) using A(r, z) obtained from our radia- 
tive transfer calculations performed on our fiducial simulation B2 (512 3 , 0.5 
Mpc). At low flux we truncate the z = 14 curve when it enters the optically 
thick domain where it no longer satisfies the ray length criterion described in 
§2.2.4. 

ization. 



which would suppress accretion of gas onto low-mass halos 
and promote homogeneity. 

The photoionization rate after reionization can be derived 
from measurements of the Lya forest. This was done by 



Kuhlen & Faucher-Giguere (2012 ) who list the values of V, 
A, and ri\ on for redshifts between 2 and 6. The quoted val- 
ues at z = 6 are T_ 12 < 0.19 and n ion < 2.6 x 10 50 s" 1 Mpc" 3 . 
Looking at Figure [7] we see that our z = 6 curve predicts an 
emissivity an order of magnitude too large for the provided 
photoionization rate. However, as we discuss in §5, our re- 
sults at z < 10 are hindered by the omission of radiative feed- 
back that would otherwise smooth inhomogeneities on small 
scales. Photoheating suppresses structure growth within the 
IGM and corresponds to a higher ionization rate at fixed emis- 
sivity due to a reduction in the recombination rate. Including 
photoheating would therefore shift the z = 6 curve in Figure 
[JJ upwards and bring us closer in agreement to the observed 
opacity of the Lya forest. At higher redshift, before photo- 
heating is important, our results should be correct, though at 
high redshift reionization becomes patchy, making a descrip- 
tion in terms of an IGM with a single UVB flux and mean free 
path less accurate. 

The strong scaling relations observed here suggest that 
small changes in n lo a can boost T by substantial amounts. Mc- 
l) use this to arg ue that the rapid evolution 
R 



2011 



Quinn etal 
in 1' observei 

by a small change in the emissivity of the ionizing background 
rather than attributing this effect to the overlap phase of reion- 



by Fan et al. ( 2006 1 at z ~ 6 can be explained 



3.5. Relationship Between Mean Free Path and Clumping 

Factor 

Suppose there is an infinitesimally thin slab of width ds 
whose area element dA is exposed to some flux F. In ion- 
ization equilibrium, the number of ionizations occurring per 
unit time balance the number of recombinations: 



dFdA = otBn e n H [[dAds, 



(14) 



where dF is the attenuation of flux passing through the slab. 
Dividing both sides of this expression by dAds, substituting 
equation ([9]), and taking n e = Hhu implies 

F(l+z)/X = a B n 2 e . (15) 

Finally, taking the clumping factor defined in equation ( flT) 
and the ionized fraction x = (n e ) / (n H ) we obtain: 



F(l+z)/X = cia B x 2 (n H ) 2 . 
For an ionized gas temperature of T gas = 10 4 K, 



A: 



23 Mpc 

X 2 C[ 



r_ 

0.3 10 



■12 „-l 
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(16) 



(17) 



where we have made use of equation (|4]) for the conversion 
from flux to photoionization rate. 

In Figure^] we plot q versus A and x 2 ci versus A from our 
radiative transfer calculations at z = 10 with T_i2 = 0.3 for 
each of the simulations listed in Table [T] In each panel the 
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FIG. 8. — c; versus A (top panel) and x 2 c/ versus A (bottom panel) for all 
simulations taken at z. = 10 with a photoionization rate r_[i = 0.3. In each 
panel points denote simulation values according to the l egen d in the top right 
corner of the plot. Traced by the dotted line is equation j l7) which describes 
the relationship between mean free path and clumping factor for an optically 
thinlGM. 

dotted line traces equation ( [T7j ) that we would expect for an 
optically thin IGM exposedto the given flux. The overall 
agreement between A and the simulation points in the bot- 
tom panel indicates consistency in the definition of the mean 
free path and detailed balance between absorptions and ion- 
izations. The minor deviations between the data points and 
the dotted curve arise because our rays have a finite length, 
s, so the mean free path evaluated in equation (TTVl will not 
correspond exactly to equation (jTOj, which is strictly correct 
only in the limit where s — > 0. 

The simulations span a large range of values in Figure [8] 
with c\ ranging from 2.6 to 8.8 and A from 2.9 to 9.5 Mpc. 
This spread arises from the broad variation in spatial and mass 
resolution exhibited by the suite of simulations. In spite of 
this, there is a clear grouping of points at c; ~ 8 and A ~ 3 
Mpc. This reflects the trend towards numerically converging 
to the "correct" clumping factor and mean free path and is our 
next topic of focus. 

4. NUMERICAL CONVERGENCE 

In this section, we attempt to answer the following ques- 
tion: What mass resolution and box sizes are necessary in 
cosmological hydrodynamics simulations, in order to obtain 
accurate results for the inhomogeneity of an unheated IGM? 
Clearly, simulations must have sufficient mass resolution to 
resolve the internal structure of the lowest mass halos that 
can contain gas. In addition, however, such simulations must 
cover a large enough volume to contain a representative sam- 
ple of the low-mass halos that dominate the opacity of the 
IGM. 

In Figure [9] we plot the clumping factors and mean free 
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FIG. 9. — Clumping factor and mean free path versus simulation box size 
(left panels) and dark matter particle mass (right panels) for z = 10 and T_i2 = 
0.3. In each panel point types distinguish between different simulations and 
follow the legend in Figure [8] Lines in the left panels connect simulations 
with fixed particle mass, while lines in the right panels connect simulations 
with fixed volume. The converged simulation C3 (1024 3 , 1 Mpc) is denoted 
by a green cross and is highlighted by a black square for easy identification. 
In the two left panels, open cyan circles and closed cyan circles denote the 
values of c; and A obtained from different random realizations of simulations 
Bl and B2 respectively. 



paths obtained from the radiative transfer calculations per- 
formed on each of the simulations listed in Table Q] To ob- 
tain a picture of convergence we display how c\ and A vary as 
functions of simulation box size, L, (left panels) and dark mat- 
ter particle mass, m& m , (right panels) for our fiducial redshift 
z = 10 and photoionization rate T_i2 = 0.3. Though the discus- 
sion below pertains explicitly to these fiducial values we have 
checked that the picture remains consistent for 10 < z < 20 
and 0.01 < T_i 2 < 10. 

Starting in the top right panel of Figure [9] we show how 
c/ changes as the mass resolution of the simulations is var- 
ied. As expected, the clumping factor increases as the reso- 
lution is refined and begins to plateau to a common value of 
c/ ~ 8.8 at mdm ~ 30 M , when a sufficiently low particle 
mass required to resolve the smallest gaseous structures in the 
box is reached. Simulations with larger particle mass are un- 
able to resolve the smallest inhomogeneities contributing to 
the dumpiness of the IGM and consequently imply clumping 
factors up to a few times smaller than the converged result. In 
the opposite limit we find that simulations with m& m < 5 M Q 
are also yielding smaller clumping factors of c; ~ 8. This ap- 
pears counterintuitive at first glance since these simulations 
should have no problem resolving the Jeans scale of the IGM. 
However, their inability to converge on q is attributed to their 
small box size, as explained below. The three smallest simu- 
lations with L = 0.25 Mpc are connected by a dot-dashed line 
and all display conspicuously unconverged values. 
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The dependence of clumping factor on box size is shown in 
the top left panel of Figure [9] In this case we expect to find 
a minimum box size above which simulations converge on q. 
Smaller boxes fail to capture power from large-scale modes 
and should therefore have reduced clumping factors. Larger 
boxes that are able to capture a representative collection of 
absorption systems within their volume should converge pro- 
vided that they have sufficient mass resolution. These general 
trends can be identified by the behaviour shown in the top 
left panel of Figure 19] There is a clear convergence of points 
near L ~ 1 Mpc andc/ ^8.8 with c\ falling off on either side. 
Simulations with L < 1 Mpc have insufficient volumes to con- 
verge with this result, while those with L > 1 Mpc have insuf- 
ficient mass resolutions. Convergence occurs in the middle 
ground where both a sufficient volume and mass resolution 
are attained. 

We find similar behaviour by comparing how A changes be- 
tween simulations. From the bottom panels of Figure [9] we 
see that the resultant behaviour is essentially an inversion of 
that described for q. Firstly, simulations with coarse mass 
resolutions overestimate A. These runs are unable to resolve 
small-scale inhomogeneities and the degree of self-shielding 
that would otherwise inhibit the propagation of ionizing pho- 
tons through neutral patches of the IGM. Secondly, the small- 
est boxes also produce values of A that are too large. As men- 
tioned above, these simulations underproduce the collection 
of halos that shield against the propagation of an ionization 
front through the IGM. 

The simulations shown here have relatively small volumes 
in a cosmological context. One issue that must be considered 
with these small boxes is that of sample variance. In the left 
panels of Figure [9] we plot c\ and A obtained from simula- 
tions Bl (512 3 , 0.25 Mpc) and B2 (512 3 , 0.5 Mpc) that were 
run using two different random realizations of the same ini- 
tial density field. The clumping factors over all three random 
realizations vary by 9% and 5% for Bl and B2 respectively, 
indicating that sample variance is somewhat important within 
these volumes. This may explain the unexpected result that 
Ci for simulation CI (1024 3 , 0.25 Mpc) is smaller than that 
of simulation Bl. Normally, at fixed volume, increasing the 
mass resolution should enhance the clumping factor (e.g., as 
seen by comparing q between simulations A3, B3, and C3). 
However, we may not expect to observe this trend with only 
one realization of a small box with large sample variance, and 
must also keep this in mind when interpreting the results of 
our convergence test. 

The above results indicate that convergence in c; and A is 
attained by simulation C3 (1024 3 , 1 Mpc) with m dm = 31 M Q . 
This simulation has a fine enough mass resolution to resolve 
small-scale inhomogeneities, and has a large enough volume 
that sample variance should be unimportant and large-scale 
modes should be captured. In order to make this claim more 
rigorous we would have to compare against a box with larger 
volume and finer particle mass than currently feasible. Never- 
theless, the data presented in Figure[9]provides compelling ev- 
idence that numerical convergence is being approached, and 
we suspect that deviations in our values for c/ and A from their 
"true" values should be small. Based on this we find that the 
necessary requirements for describing the inhomogeneity of 
an unheated IGM using cosmological hydrodynamics simula- 



tions is to use box sizes L > 1 Mpc with dark matter particle 
masses m<i m < 50 Mq. Smaller boxes are troubled by sample 
variance while coarser mass resolutions are unable to resolve 
the mass scale where gaseous halos are dominating the opac- 
ity of the IGM. 

5. DISCUSSION 

We have performed high-resolution, cosmological simula- 
tions of structure formation at redshifts z > 6, including adia- 
batic hydrodynamics. By post-processing the resulting den- 
sity fields with a radiative transfer algorithm for hydrogen 
ionizing radiation, we have determined the opacity of the un- 
heated IGM, in terms of the mean free path to ionizing ra- 
diation, A, as a function of redshift and ionizing background 
intensity. These results are relevant (1) as converged solutions 
for the opacity of the IGM early in the reionization process, 
before photoheating has evaporated small-scale structure and 
(2) in determining what mass and length resolutions are nec- 
essary to correctly model the propagation of ionization fronts 
into the neutral IGM. We derive values of n cl ; t , the proper hy- 
drogen number density above which gas remains neutral, that 
are for the most part a function of only T_i2. Simulations that 
mimic the effect of self-shielding by turning off the optically- 
thin flux at high densities should use n cr ; t ^0.1 cm" 3 TtQ, 
independent of redshift. 

Our post-processing approach neglects the hydrodynamic 
feedback of photoheating on the density evolution. These re- 
sults therefore indicate what the initial degree of inhomogene- 
ity should be as ionization fronts propagate into the IGM. In 
addition, they place an upper limit to this inhomogeneity in 
patches of the IGM that have already been ionized. We find 
that the initial clumping factor of the IGM just as it is be- 
ing ionized is a strong function of redshift and ionizing back- 
ground intensity, with typical values at z = 10 ranging from 
about q = 4.4 to 16 and A = 0.7 to 15 Mpc, for r_i 2 = 0.03 to 
r_i2 = 3, respectively. 

Modelling the transition from a neutral to ionized IGM 
requires self-consistent simulations of the coupled radia- 
tive transfer and hyd rodynamical photoevaporation process. 
Shapiro et al. (2004) used idealized two-dimensional radia- 
tive transfer hydrodynamics calculations of the photoevapo- 
ration of initially spherical, isolated minihalos, surrounded by 
infalling gas. Those calculations showed that smaller mini- 
halos are photevaporated faster, and that larger fluxes lead to 
faster photoevaporation times as well. However, the photoe- 
vaporative process is in reality likely to be more complex than 
for the simplified geometries and source lifetimes considered 
by |Shapiro et al.| ( |2004| l, with halos over a range of masses 
clustered in space and arranged within a "cosmic web" of fila- 
mentary structure. Filamentary infall from nearby neutral gas 
could replenish halos as they are being evaporated, consider- 
ably extending the photoevaporation process, while ionization 
from highly luminous but intermittent starburst galaxies could 
result in large clumping factors and stalled minihalo evapora- 
tion, considerably increasing photon consumption and leading 
to a much more com plex morphology of early H II regions 
(e.g., Wise et al.|2012| l than is typically envisioned. 

We find that convergence is reached at a dark matter particle 
mass of mdm < 50 Mq. A box size of L > 1 Mpc is necessary 
to sample the IGM for the purpose of modelling absorptions 
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FIG. 10. — Here we compare clumping fa ctors C-i fr om our fiducial sim- 
ulation B2 (512 3 , 0.5 Mpc) to the |Pawlik etld"|j2009) simulation L6N256 
(256 3 , 6 Mpc). The Pawlik et al. curve corresponds to their reference sim- 
ulation that did not include photoionization and has similar attributes to our 
simulation A6 (256 3 , 8 Mpc). As expected, A6 and L6N256 are in good 
agreement with each other while the higher-resolved simulation B2 shows 
clumping factors a few times larger than L6N256. 



by small-scale structure. The clumping factors we find from 
our converged results are somewhat smaller than the values 
q ~ 30 found in early attempts to characterize the dumpiness 
of the IGM whi ch did not accurately sepa rate ionized and neu- 
tral gas (e.g., |Gnedin & Ostrikerf l997), but are higher than 
the clumping factors found by Pawlik et al. (2009) at z ~ 9, 
just before the IGM in those simulations was heated by ion- 
izing radiation. We attribute this difference to the increased 
mass resolution of our simulations, which resolve halo masses 
down to the Jeans mass in an untreated IGM (~ 10 4 M Q ), as 
opposed to that corresponding to the Jeans mass for a pho- 
toioniz ed gas temperature o f ~ 10 4 K (~ 10 8 M Q ). As pointed 
out by Pawlik et al. (2009), the clumping factors they find at 
z ~ 6, for a patch or the IGM which was ionized significantly 
earlier, at z ~ 9, are likely to be correct, because a long enough 



time had passed since the gas was ionized for photoheating to 
evaporate existing small-scale structure and suppress accre- 
tion onto newly formed dark matter minihalos with masses 
below - 10 8 " 9 M Q , which were resolved in their highest res- 
olution simulations by ~ 100- 1000 dark matter particles. 

Here we point out that their values of the clumping fac- 
tor are likely underestimates before and for sometime after 
reionization, by about a factor of a few. This is illustrated 
in Figure 10 in which we show our clumping factors, C_i, in 
an unhe ated IGM for n^ t = 0.1 cm -3 , along with those found 
by Pawlik et al. (2009i for the same threshold oversensitM 3 ] 
The clumping factors we find in our fiducial simulation B2 
are at all times larger, by a factor of 1.2 at z = 20 and 3.5 at 
z = 6. As described above, the clumping factors we obtain at 
low redshift are overestimates, and in reality should be closer 
to C-i ~ 6 found by Pawlik et al. in the presence of a pho- 
toevaporative background. Combining these two results, we 
find that the clumping factor of the IGM evolves strongly just 
after a patch of IGM is ionized. For ionization at z = 9, the 
clumping factor drops from c\ ~ 20 at z = 9 to a few at z = 6, 
depending on the intensity of the ionizing background - with 
larger intensities leading to higher clumping factors and larger 
mean free paths. The Pawlik et al. clumping factor shown 
in Figure [TO] corresponds to their unheated simulation named 
L6N256 (256 3 , 6.25 Mpc), which is similar to our simulation 
A6 (256 3 , 8 Mpc). Indeed, our simulation A6 does a good job 
of reproducing their result. 

The results presented here for the inhomogeneity of elec- 
tron density in the presence of an ionizing background should 
serve as a foundation for more detailed study of radiative 
transfer and hydrodynamical effects in the initial stages of 
reionization, including the effects of the initi al relative veloc- 
ity between baryons and dark matter (e.g., Tseliakho vich &| 
|Hirata i2010), pr eheating by long mean free path X-ray pho- 
tons (e.g., |Ricotfi & Ostriker||2004[ |Ricotti et al.||2005], and 
photoevaporation (e.g., Shapiro et aL|2004||Abel et al.|2007[ ). 
In simulating all of these processes, it will be necessary to 
resolve small-scale structure in the way outlined here. 

We thank P. R. Shapiro for helpful discussions. The simu- 
lations used in this paper were performed on the GPC super- 
computer at the SciNet HPC Consortium. SciNet is funded 
by: the Canada Foundation for Innovation under the auspices 
of Compute Canada; the Government of Ontario; Ontario Re- 
search Fund - Research Excellence; and the University of 
Toronto. 
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